#include <stdlib.h>
#include <fstream>
#include "Defines.h"

main()
{
    printf("Test\n");

    Solver pSolver; // the solver object
    
    char solvername; // name of the solver to use
    int Nx1, Nx2; // dimensions of the domain
    REAL psitol; // tolerance of the solver
    
    Array2D<REAL> Sc(); // center stencil
    Array2D<REAL> Ss(); // south stencil
    Array2D<REAL> Sw(); // west stencil
    
    Array2D<REAL> rhs(); // right hand side
    Array2D<REAL> sol(); // exact solution
    Array2D<REAL> est(); // initial guess of solution

    int ii,jj; // counting variables
    
    // Define the type of solver
    solvername = 'rrb';
    
    // Define dimension of the problem
    Nx1 = 200;
    Nx2 = 200;
    
    // Define tolerance of the iterative solver
    psitol = 2e-6;

    // Define the center, south and west stencil
    Sc.resize(Nx1+3, Nx2+3);
    Ss.resize(Nx1+3, Nx2+3);
    Sw.resize(Nx1+3, Nx2+3);

    Sc = 0.0;
    Ss = 0.0;
    Sw = 0.0;
    
    // Define the right hand side and initial solution
    rhs.resize(Nx1+3, Nx2+3); // right hand side
    sol.resize(Nx1+3, Nx2+3); // exact solution
    est.resize(Nx1+3, Nx2+3); // initial guess of solution

    rhs = 0.0;
    sol = 0.0;
    est = 0.0;

    for (jj = 1; jj <= m_Nx2; ++jj)
    {
        for (ii = 1; ii <= m_Nx1; ++ii)
        {
            sol[ii][jj] = 1.0 + 10*((REAL)(ii + jj))/(Nx1 + Nx2);
        }
    }

    for (jj = 1; jj <= m_Nx2; ++jj)
    {
        for (ii = 1; ii <= m_Nx1; ++ii)
        {
            rhs[ii][jj] =  Ss[ii    ][jj + 1] * sol[ii    ][jj + 1]
                         + Sw[ii + 1][jj    ] * sol[ii + 1][jj    ]
                         + Ss[ii    ][jj    ] * sol[ii    ][jj - 1]
                         + Sw[ii    ][jj    ] * sol[ii - 1][jj    ]
                         + Sc[ii    ][jj    ] * sol[ii    ][jj    ];
        }
    }
    
    // Create a solver
    pSolver = SolverFactory::create(solvername, &Sc, &Ss, &Sw, Nx1, Nx2, psitol);

    // Call the solver
    pSolver -> solve(&rhs, &est);


}

